1 Preparation

1.1 Load packages

library(sf)          # classes and functions for vector data
## Warning: package 'sf' was built under R version 3.5.3
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(sp)
library(spData)        # load geographic data
## Warning: package 'spData' was built under R version 3.5.3
library(spDataLarge)   # load larger geographic data
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(stplanr)      # geographic transport data package
## Warning: package 'stplanr' was built under R version 3.5.3
library(mapview)
## Warning: package 'mapview' was built under R version 3.5.3
library(readxl)
library(xlsx)
## Warning: package 'xlsx' was built under R version 3.5.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.3

1.2 Set up API keys

Sys.setenv(CYCLESTREET = "4f8888f9c99f438b")
Sys.setenv(GRAPHHOPPER='4897ffda-408e-4f6b-8161-21ecce6ae1d1')

2 Import Data

ct_flow<- read_excel("GEOG506/ct_flow_island.xlsx") #commuting flow

load("Data/origin_ct.Rdata") #origin cT (residential area)
load("Data/dest_ct.Rdata") #destination CT (industrial, commercial, institutional area)
load(file="Data/CT.Rdata") #CT Island of Montreal
head(ct_flow)
head(Island)

3 Flow

zones_intra = filter(ct_flow, Origin == Destination)
zones_inter = filter(ct_flow, Origin != Destination)
flow_lines = od2line(zones_inter, origin, dest)
flow_lines
## class       : SpatialLinesDataFrame 
## features    : 62762 
## extent      : 269269.1, 305880.3, 5029535, 5061994  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## variables   : 7
## names       :     Origin, Destination, Total, Car, Transit, Walking, Cycling 
## min values  : 4620001.00,  4620002.00,     0,   0,       0,       0,       0 
## max values  : 4620619.00,  4620619.00,   975, 330,     485,     440,      45

4 Generate Routes

flow_gcs <- st_transform(st_as_sf(flow_lines), "+proj=longlat +datum=WGS84") #transform coordinate system
#generate car routes using graphhopper
#generate bicycle routes using cyclestreet
Car = line2route(flow_gcs, "route_graphhopper")
cycling = line2route(flow_gcs, route_cyclestreet)
Car
cycling

5 Emission Saving

5.1 Potentially Cyclable Trip

cycling_routes <- as(cycling,"Spatial")
car_trips = left_join(Car,cycling_routes@data[,c("time","length","id","up_tot")], by = "id")
car_trips
#calculate time ratio
car_trips$ratio = car_trips$time.y/60/car_trips$time.x
car_trips
#scenario 1：cycling distance <=4.6 km, elevation gain <= 30 m
car_1 <- car_trips[car_trips$length <= 4600 & car_trips$up_tot <= 30,]
#scenario 2：cycling distance <=4.6 km, elevation gain <= 30 m, ratio <=1.2
car_2 <- car_trips[car_trips$length <= 4600 & car_trips$up_tot <= 30 & car_trips$ratio <=1.2,]

5.2 Emission Saving by Origin

5.2.1 Scenario 1

flow_lines@data$id <- row.names(flow_lines@data)
car_1 <- left_join(car_1,flow_lines@data[,c("Origin","id","Car")], by = "id")
#calculate emission saving per flow
car_1$Saving <- (car_1$dist/1000000)*car_1$Car*19
car_1
save_1 <-car_1 %>% 
  group_by(Origin) %>% 
  summarize_if(is.numeric, sum) %>% 
  dplyr::rename(CTUID = Origin)
save_1_sp <- as(save_1,"Spatial")
save_1_joined = left_join(Island,save_1_sp@data[,c("CTUID","Saving")], by = "CTUID")

5.2.2 Scenario 2

car_2 <- left_join(car_2,flow_lines@data[,c("Origin","id","Car")], by = "id")
#calculate emission saving per flow
car_2$Saving <- (car_2$dist/1000000)*car_2$Car*19
car_2
save_2 <-car_2 %>% 
  group_by(Origin) %>% 
  summarize_if(is.numeric, sum) %>% 
  dplyr::rename(CTUID = Origin)
save_2_sp <- as(save_2,"Spatial")
save_2_joined = left_join(Island,save_2_sp@data[,c("CTUID","Saving")], by = "CTUID")

6 Maps

6.1 Cyclable Route

6.1.1 Scenario 1

car1_ln = SpatialLinesNetwork(car_1)
edge_best1 = igraph::edge_betweenness(car1_ln@g)
mv1 = mapview(car1_ln@sl$geometry,lwd=edge_best1 / 3000,color = "chocolate", alpha = 0.8)
mv1@map

Cyclable Trips (Scenario 1)

6.1.2 Scenario 2

car2_ln = SpatialLinesNetwork(car_2)
edge_best2 = igraph::edge_betweenness(car2_ln@g)
mv2 = mapview(car2_ln@sl$geometry,lwd=edge_best2 / 3000,color = "purple", alpha = 0.8)
mv2@map

Cyclable Trips (Scenario 2)

6.2 Emission Saving

6.2.1 Scenario 1

mv3 = mapview(save_1_joined,zcol = ("Saving"),alpha = 0.8, legend = TRUE, layer.name="Emission Savings (kg)")
mv3@map

Emission Savings Per Day (Scenario 1)

6.2.2 Scenario 2

mv4 = mapview(save_2_joined,zcol = ("Saving"),alpha = 0.8, legend = TRUE, layer.name="Emission Savings (kg)")
mv4@map

Emission Savings Per Day (Scenario 2)